#!/usr/bin/env python
# coding: utf-8

# if in rstudio: reticulate::use_condaenv('scvelo-colorbar')

# # Getting Started
# input data for scVelo are two count matrices of pre-mature (unspliced) and mature (spliced)
# abundances, which can be obtained from standard sequencing protocols, 
# using the `velocyto` or `kallisto` counting pipeline.
# 
# ## workflow at a glance 

# import scvelo
import scvelo as scv

# change matplotlib setting to scv's default
scv.set_figure_params(style='scvelo')

# * default is `scvelo` style, no border or grid in plot
# 
# * set to `scanpy` style to add border and grid to plot
# 
# * set to `None` is very similar to `scanpy`?

# ## Read your data

# can read h5ad, loom, csv
# set cache=True to speed up reading by creating cache file on disk
adata = scv.read('/home/gjsx/learn/scvelo/data/Pancreas/endocrinogenesis_day15.h5ad', cache=True)

# can merge spliced & unspliced data
# scv.utils.merge()

# ### expression matrix
adata.X

# ### annotation of cells / observations
adata.obs


# ### annotation of genes / vars
adata.var


# ### unstructured data like umap
adata.uns


# ### data layers where spliced and unspliced counts are stored
adata.layers


# ## Basic preprocessing

# basic preprocessing (gene selection and normalization)
# keep only genes with counts>=20 in spliced+unspliced data
# select high variable 2000 genes for pca
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
adata

adata.obs

# compute the first- and second-order moments (means and uncentered variances)
# for velocity estimation
scv.pp.moments(adata, n_pcs=30)
# Ms and Mu layers are added (Moment-spliced + Moment-unspliced)
adata

# This is the core function of scvelo
# default mode is stochastic (Bergen2019)
scv.tl.velocity(adata, mode='stochastic')
# dynamic model (Bergen2019) need scv.tl.recover_dynamics(adata) beforehand.

# velocity are stored in layers
adata.layers

# deterministic modal (Monne2018) also can be used
adata_deter = scv.tl.velocity(adata, mode="deterministic", copy=True)
# deterministic model calc no variance of velocity
adata_deter

# ### project velocity to low-dimension

scv.tl.velocity_graph(adata)

# * A sparse matrix of velocity graph (cellnumber\*cellnumber) store to uns
# * similarities between every cell's timely states?

adata.uns


# ## Visualization
# ### draw vector as small points (not very pretty)

scv.pl.velocity_embedding(adata, basis='umap')


# ### draw stream-like long arrow (with cell type texts)

# change color param to change annotation
scv.pl.velocity_embedding_stream(adata, basis='umap', color="clusters")


help(scv.pl.velocity_embedding_stream)


# ### draw velocity arrow in grid

scv.pl.velocity_embedding_grid(adata, basis='umap')

# ### counterparts of tools module (scv.tl)

# examine specific genes velocity
scv.pl.velocity(adata, var_names=['Cpe',  'Gnao1', 'Ins2', 'Adk'], ncols=2)


# * every gene has a velocity value in a cell.
# * positive means upregulating, negative means downregulating
# * then the dim-reduc projection is just use velocity value instead of expression level?

scv.pl.velocity_graph(adata)

